{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Computing the SVD"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "np.random.seed(15)\n",
        "n = 5\n",
        "A = np.random.randn(n, n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now compute the eigenvalues and eigenvectors of $A^TA$ as `eigvals` and `eigvecs` using `la.eig` or `la.eigh` (symmetric):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "eigvals, eigvecs = la.eigh(A.T @ A)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([0.08637178, 0.457892  , 2.04177547, 2.34383161, 8.37000184])"
            ]
          },
          "execution_count": 12,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "eigvals"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Eigenvalues are real and positive. Coincidence?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(5, 5)"
            ]
          },
          "execution_count": 13,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "eigvecs.shape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Check that those are in fact eigenvectors and eigenvalues:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1.77635684e-15,  2.22044605e-15,  3.10862447e-15,\n",
              "         8.88178420e-16,  8.88178420e-16],\n",
              "       [ 2.22044605e-16, -8.88178420e-16,  0.00000000e+00,\n",
              "        -8.88178420e-16,  2.63677968e-16],\n",
              "       [-8.88178420e-16,  4.44089210e-16,  8.88178420e-16,\n",
              "        -2.22044605e-16,  6.66133815e-16],\n",
              "       [ 0.00000000e+00, -4.44089210e-16,  8.88178420e-16,\n",
              "         0.00000000e+00,  7.77156117e-16],\n",
              "       [ 0.00000000e+00, -9.02056208e-17,  3.88578059e-16,\n",
              "         0.00000000e+00,  3.33066907e-16]])"
            ]
          },
          "execution_count": 14,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "B = A.T @ A\n",
        "B - eigvecs @ np.diag(eigvals) @ la.inv(eigvecs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`eigvecs` are orthonormal! (Why?)\n",
        "\n",
        "Check:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "2.094158390904296e-15"
            ]
          },
          "execution_count": 15,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(eigvecs.T @ eigvecs  - np.eye(n))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "------\n",
        "Now piece together the SVD:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Sigma = np.diag(np.sqrt(eigvals))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "V = eigvecs"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "U = A @ V @ la.inv(Sigma)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Check orthogonality of `U`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[-1.11022302e-15, -7.67324533e-16, -6.67804308e-16,\n",
              "         1.47963782e-15,  1.44717415e-16],\n",
              "       [-7.67324533e-16,  6.66133815e-16, -5.43672425e-16,\n",
              "         6.51481016e-16, -2.64998805e-16],\n",
              "       [-6.67804308e-16, -5.43672425e-16,  2.66453526e-15,\n",
              "        -1.42306617e-15,  1.55464267e-15],\n",
              "       [ 1.47963782e-15,  6.51481016e-16, -1.42306617e-15,\n",
              "         3.10862447e-15, -1.13201575e-15],\n",
              "       [ 1.44717415e-16, -2.64998805e-16,  1.55464267e-15,\n",
              "        -1.13201575e-15,  4.44089210e-16]])"
            ]
          },
          "execution_count": 19,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U @ U.T - np.eye(n)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.2+"
    },
    "widgets": {
      "state": {},
      "version": "1.1.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}